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SUMMARY 

We adopt the frozen Gaussian approximation (FGA) for mod- 
eling seismic waves. The method belongs to the category of 
ray-based beam methods. It decomposes seismic wavefield 
into a set of Gaussian functions and propagates these Gaussian 
functions along appropriate ray paths. As opposed to the clas- 
sic Gaussian-beam method, FGA keeps the Gaussians frozen 
(at a fixed width) during the propagation process and adjusts 
their amplitudes to produce an accurate approximation after 
summation. We perform the initial decomposition of seismic 
data using a fast version of the Fourier-Bros-Iagolnitzer (FBI) 
transform and propagate the frozen Gaussian beams numeri- 
cally using ray tracing. A test using a smoothed Marmousi 
model confirms the validity of FGA for accurate modeling of 
seismic wavefields. 



INTRODUCTION 

Ray t heory ( ICervenvll200ll : |PoDovll2002luEngquist and Runbord . 
120031) is a widely used approach to seismic modeling and mi- 
gration. In this approach, one decomposes the wavefields into 
elementally waveforms that propagate along rays, and then con- 
structs Green's functions or wavefields according to the dy- 
namic information on rays (e .g. path traje ctory, amplit ude and 
phase ). Kirchhoff migration JGravl 1 1 98a ; j Keho and B evdounl. 
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1988") and Gaussian beam migration ('Hill,'l990, 2001 ; Nowack et 
2003; Grav, 2005; Gray and Bleistein, 2009; Popov et al., 2010) 
are the most famous seismic applications of this principle. 

Dynamic ray tracing used in Kirchhoff migration is an effec- 
tive method, however it produces unbounded amplitudes at 
caustics. The Gaussian beam approximation (GBA) retains the 
merits of ray tracing but can also handle multipathing while 
maintaining accuracy at caustics. However, in order to obtain 
a good resolution in GBA, one needs to tune the width param- 
eter of Gaussian beams, especially because be ams spread sig- 
nificantly in the process of wave propagation JCervenv et al.L 
ll982l ; lHnl ll990'; Fomel an d Tanushevl 2009). This parameter 
tuning becomes difficult in practical applications because of 
the heterogeneity of the media and the non-linearity of the Ric- 
cati equation involved in the beam construction. GBA relies on 
a Taylor expansion around the central ray, hence the error of 
the approximation in creas es when the beams become wider. 
lOian and Yiii3 ( I2OIOI) and ILu and Yand ( 1201 ih analyzed this 
phenomenon and showed that, even in the case of a simple ve- 
locity distribution, the error of GBA may grow rapidly in time. 

In this paper, we adopt the frozen Gaussian approximation 
(FGA) method for computing seismic wave propagation in com- 
plex structures. The method was introduced in previo us stud- 
ies on genera l linear strictly hyperbolic systems (.Lu and Yang . 
l201lLl2012al lbl) and was originally motivated by Hetman-Kluk 



propagator for solving the Schrodinger equation in qu antum 
chemistrv (" Herman and Klukl . I I984I ; Kavl . 1 1 9941 . [200a) . The 
main idea of FGA is to use Gaussian functions with fixed widths 
to approximate the propagation of seismic wavefields. These 
Gaussian functions are also known under the name coherent 
states in quantum mechanics, and mathematically form a tight 
frame in L -function space. The coherent state method was 
previously a pplied in seismic imaging (Albertin et al., 2001; 
[Foster et al.l 12002) but lacked the rigorous treatment of ampli- 
tude factors provided by FGA. Despite its superficial similarity 
with GBA, FGA is different at a fundamental level: In GBA, 
on the other hand, each Gaussian beam provides an approxi- 
mate solution to the wave equation, which requires its width to 
change in time l IBleistein and GravluOlOr) . In FGA, Gaussian 
functions are used only as building blocks for wave propaga- 
tion. Each individual frozen Gaussian does not approximate a 
solution to the wave equation. Compared with GBA, FGA is 
able to provide a more accurate and robust solution due to the 
error cancellatio n phenomenon, especially in the situation of 
wave spreading dLu and Yand.l201Ul2012ah . 

The main procedure of FGA is the following. The initial data 
are decomposed into a set of Gaussians with fixed predefined 
naiTow widths. Each Gaussian function is then propagated 
along geometric rays. The amplitudes of the Gaussians are ad- 
justed according to rigorously-derived dynamic equations so 
that, at the final time, the sum of them yields an accurate ap- 
— .proximation to seismic wavefields or Green's functions. The 
^ 'basic ingredients of the FGA implementation are similar to 
those used in other ray-type or beam methods. 

To accelerate numerical computation, we develop an algorithm 
for a fast Fourier-Bros-Iagolnit zer (FBI) transform algorithm 
l IFollandll 1989l ; lMartinezll2002l) to decompose the initial wave- 
fields into Gaussians. The idea of decomposing functions into 
Gaussians is k nown i n signal pro cessing applications as Ga- 
bor expansion (lGaboil.f r946, ; .Helstromlll966l ; lBastiaansl[T980l ; 
iMa and Marg rave. 2008). Thresholding criteria of initial Gaus- 
sian packets can be chosen to provide different levels of details 
depending on accuracy requirements and computational cost 
constraints. 



THEORY 

For simplicity, we shall consider wave propagation governed 
by the acoustic wave equation in d dimensions. The seismic 
wavefield u^ is a function of time t and spatial variable x e M , 



with initial conditions 

u'{Q,x)=f^{x), 
d,u'(Q,x)=fl[x), 



(1) 

(2) 
(3) 
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where c(x) is the seismic velocity, and e indicates the depen- 
dence on the (rescaled) wave length. We assume that e <S 1 is 
the high-frequency wave regime corresponding to short wave- 
lengths. 

FGA Formulation 

FGA is an asymptotic-based approach to high-frequency wave 
propagation. Its complete derivation, error estimates (includ- 
ing validity at caustics), and generalization to other strictly hy- 
perbolic systems are provided bv Lu and Yang (2011, 2012a .b). 

FGA gives the asymptotic approximation to the solution of 
equation lUll as a sum of Gaussian functions with fixed width. 
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where 

Viiq,p)= I ul^,{y,q,p)e--My-^)-T.\y-^\' Ay, 

4,o(v,.,P)4(/,fW±^/fW 



(5) 
(6) 



In equation ([4}, i = %/— T is the imaginary unit, and "+" and 
"— " indicate the two wave branches, and G± are th e sets of 
(g, p) pairs. Equation (|5} is in a form of FBI transform ( iMartinea. 
|2002). In FGA, associated with each frozen Gaussian, the 
time-dependent quantities are: the center 2^, momentum P± 
and amplitude a±. The weight function \ff± is time-indepen- 
dent and computed initially, while the width of the Gaussian 
function is fixed at all times. 

The evolution of Q^(t,q,p) and P±{t,q,p) satisfies the ray 
tracing equations corresponding to the Hamiltonian 



H± = ±ciQ±)\P±\ 



(7) 



For ease of notations, we will suppress the subscripts "±" 
when no confusion might occur. Hence (Q,P) solves 



dg 

dt 
dP 



dpH, 



-SqH, 



with initial conditions 

Q{0,q,p)=q and P{Q,q,p)=p. 

The evolution of the amplitude a{t,q,p) is given by 

da dpHdnH a ( , dZ\ 

— =a— ^— + -tr Z"' — . 

df H 2 \ dt J 



(8) 



(9) 



(10) 



with initial condition a{0,q,p) = 2''". In dlOt , we have used 
the shorthand notations 



d, = dg - id, 



P' 



Z = d-SQ + iP). 



(11) 



Here d-Q and d^P are understood as matrices, with the {j.k) 
component of a matrix d-Q given by d,.Qi^. The matrices d-Q 



and d,P can be solved at each time step by either divided dif- 
ference or the following dynamic ray tracing equations. 



did,P) ^ d^H ^ d^H 

= -d-,Q ^-d.P^-^:—. 

dt dQ^ ~ 3PdQ 



(12) 



(13) 



Componentwise, (I12l(-J13l> can be written as (with Einstein's 
index summation convention) 

z^ijK ^d^ Q — .^ ^d^ Pi 



dt 

d{d:P)jk 
dt 



-dz,Qi 



dQidPk 



dQiQk 



d,^P, 



dPiPk' 
dPldQk 



We need to point out a key difference here: In FGA, the solu- 
tion of dynamic ray tracing equations ( I12ll3t only affects the 
amplitude a, while in GBA, it affects both the amplitude and 
the beam width. We also note that equation dlOl ) actually gives 
a conserved quantity along the Hamiltonian flow l[8j. 



c2(e)detZ 



= constant. 



which implies that 



a='^id..Zr\ 

c{q) 



(14) 



(15) 



with the appropriate branch of the square root continuously 
determined in time. 

Initial wavefield decomposition 

In order to apply representation (|4) in practice, it is necessary 
to decompose the initial wavefield l|2j-l|3j into a sum of Gaus- 
sian functions, i.e. , to choose proper sets G± of {q,p) pairs in 
equation Q and compute t//ji in equation ^ correspondingly. 
Here we propose a local fast FBI transform to efficiently com- 
pute l//^. 

Let us rewrite equation (|5} in the form 

Vj{q,p)= I f'j{y)e-'iP''^y-''^-TAy-i\'dy, (16) 

for J = 0, 1. Equivalently, 

Vj{q,P)= f f^{q + r)e-iP-'-TMdr, (17) 

where we use the change of variable r = y — q. Define 



8lj{r)=.ff{q + r)exp{-^^\r\') 



(18) 



then t//? is given by the (rescaled) Fourier transform of g^ ■, 

Vj{g,p)=glj{p/£)- (19) 

Notice that g^ contains an exponential function, hence its 
function value is negligible outside a localized domain cen- 
tered around zero, for example, a small box 

Se = [-L/2,L/2]''cR'^ 
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with the length L scaled as ff{-\/e). Therefore, \I/J{q,p) can 
be evaluated efficiently by applying Fast Fourier Transform of 
g^ . restricted on the small box B^. 

Once y/jL is computed, one can apply a simple thresholding to 
select the sets G± where y/^ have relatively large values. 

Algorithm 

In summary, the FGA algorithm consists of three steps: 

1. Initial decomposition: Choose the sets G± of {q,p) 
pair and calculate t/Zji defined in equation ^ from /g 
and/f; 

2. Time propagation: Numerically integrate (|8) and JIO) 
up to the final time T; 

3. Reconstruction: Compute the wave field at time T by 
applying equation l|4ll. 



With the modifications proposed bv lTanushev et alj ( 1201 11) . it 
is possible to generate initial data for Step 1 from the seis- 
mic wavefield recorded at the Earth surface, and to incorporate 
FGA into a seis mic imaging schem e in the fashion of reverse- 
time migration dEtgen et al.Ll2009t) . In the next section, we test 
the accuracy of the approximation itself by modeling a Green's 
function (wave propagation from a point source) in a synthetic 
velocity model. 



NUMERICAL TEST 

We test the performanc e of FGA using a smoothed Marmousi 
model jVersteegl 1 19940 shown in Figure [T] Our goal is to ex- 
trapolate the expanding wavefield from a point source from 
the initial condition shown in Figure|2]through 0.25 seconds in 
time. The reference calculation for comparison is performed 
with the high l y accu rate lowrank symbol approximation method 
' r5ll2013l) . 
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For the initial beam decomposition, we set the value of e by 
the ratio of £ norms of the initial wavefield /q and the wave 
speed f^. We adopt a fourth-order Runge-Kutta method as the 
numerical integrator for time propagation. 

The results in Figure|3]are produced by A' = 48, 521, and 5650 
frozen Gaussian beams, which corresponds to keeping {q,p) 
pairs with W-largest amplitudes of l|/^. The time step is taken 
as 0.001 s. In the initial wavefield decomposition, the mesh 
sizes for q and p are 0.08 km x 0.06 km and 0.0898 km^' x 
0.0898 km^ ' respectively, and the mesh size for y is 0.004 km x 
0.004 km. The images are reconstructed using the mesh size 
0.004 km X 0.004 km. The result in Figure|4]uses only 6 beams 
for the extrapolation. It is easy to see that the individual beam 
width remains "frozen" and does not change with position. 

Comparing the results in Figure |3] with the reference wave- 
field shown in Figure [5(a)] we can see that propagating only 
N = 521 Gaussians is sufficient to produce a qualitatively ac- 
curate result, while N = 5650 Gaussians produce a quantita- 
tively accurate result by adjusting amplitudes and generating 



small-scale wavefield features including caustics. The differ- 
ence between the result in Figure [3(c)] and the reference wave- 
field plot is plotted in Figure [5(b)] at the same scale and appears 
to have negligible magnitude. 
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Marmousi IVIodel: Smoothed Velocity 

Figure 1: Smoothed Marmousi velocity model. 
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Figure 2: Initial seismic wavefield generated by a point source 
in the center of the Marmousi model. 



CONCLUSIONS 

The Frozen Gaussian Approximation (FGA) provides a stable 
and efficient computational tool for seismic wavefield propa- 
gation. FGA uses Gaussian functions with fixed widths, rather 
than those that might spread over time, as in the classic Gaus- 
sian beam approximation. As a result, a stable behavior and 
a good approximation accuracy can be achieved without tun- 
ing the width parameters. A rigorous mathematical analysis 
guarantees the accuracy of the FGA solution in modeling wave 
propagation beyond caustics. We have also introduced a local 
fast FBI transform algorithm that decomposes the initial wave- 
field into a set of Gaussian functions. We have numerically 
tested the performance of FGA using a smoothed Marmousi 
model. The method may find direct applications in seismic 
modeling and seismic imaging by beam migration. 
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Figure 3: Wavefield predicted by FGA using (a) N = 48, (b) 
A^ = 521, and (c) A' — 5650 Gaussians for each wave branch. 
The accuracy of the approximation gets improved by adding 
more beams. 




Figure 4: Wavefield at 0.25 s after the initial wavefield from 
Figure|2l as predicted by FGA using only six initial Gaussians 
for each wave branch (A' = 6 in equation|4j. 
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Figure 5: (a) Reference wavefield. (a) Quantitative difference 
between with the result in Figure [3(c)] 
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